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Abstract 

Forty-five years after the point de depart [l] of density functional theory, its ap- 
plications in chemistry and the study of electronic structures keep steadily growing. 
However, the precise form of the "divine" energy functional in terms of the electron 
density [2] still eludes us — and possibly will do so forever j.'jj . In what follows we 
examine a formulation in the same spirit with phase space variables. The validity of 
Hohenberg-Kohn-Levy-type theorems on phase space is recalled. We study the repre- 
sentability problem for reduced Wigner functions, and proceed to analyze properties of 
the new functional. Along the way, new results on states in the phase space formalism 
of quantum mechanics are established. Natural Wigner orbital theory is developed 
in depth, with the final aim of constructing accurate correlation-exchange functionals 
on phase space. A new proof of the overbinding property of the Miiller functional is 
given. This exact theory supplies its home at long last to that illustrious ancestor, the 
Thomas-Fermi model. 

To Jens Peder Dahl, dear friend and Wigner function's stalwart 
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1 Introduction 

1.1 Conventions and perspective 

In this article Hartree atomic units [4] are used. The operator Hamiltonian EI for N fermions 
involves only one- and two-body observables. We work under the Born-Oppenheimer regime 



2 



and look exclusively at the electronic problem, regarding the potential due to presence of 
nuclei as an external one. To fix ideas, consider the problem of N electrons in an ion 
of charge Z in the common approximation that neglects spin-orbit interaction and weaker 
couplings; so that in configuration space 
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Remember that the set of all iV-particle density matrices Dn coincides with the set T> of 
positive hermitian operators of unit trace on the Hilbert space of antisymmetric iV-particle 
functions. This is a convex set, and its extreme elements are the pure states. When the 
system is in the (normalized) pure state \^n) the iV-particle density matrix is of the form 
D N = \^f N )(^f N \; then D 2 N = D N . Given D N , one refers to the integral operator with kernel 
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as the reduced n-matrix, or simply as the n-matrix. We denote by the same letter a n-matrix 
and its kernel and employ the standard notation 1 = (q^Si), and so on, for the spatial and 
spin variables]^] These reduced matrices are still positive operators with trace 1, but now 
they form only a proper subset of T> n . The case n = 2 is of special importance, since the 
energy of system takes the form 



E 



tr(KD 2 ). 



(2) 



This is an exact linear functional of D 2 , and the ground-state energy would be obtained 
by minimizing it on the set of 2-matrices. Therein lies the rub, as the iV-representability 
problem for 2-matrices has never been efficiently solved [5). This is why the Hohenberg- 
Kohn-Sham density functional theory (DFT) program and its generalizations, initiated by 
the proof that the ground state density determines every property of an electronic system, 
still enjoy a tremendous success. 



1.2 Purpose and plan of the article 

This paper has the aim of sprucing up an untended corner of the garden of generalized density 
functional theories. The original DFT program is beset by our foreordained ignorance of 
the exact energy functional, while the reduced density matrix approach is hobbled by the 
unsolvability of the iV-representability problem for two-electron densities]^] We push forward 

1 In point of rigour, the partial diagonals in ([T]) are not generally defined; however, one can always make 
sense of the formula by means of the spectral theorem. 

2 Even before the results in [3] the problem has been regarded as well nigh intractable. Nevertheless, there 
are the formal solution by Garrod and Percus [6], still an appropriate reference for TV-representability; Ayers' 
reformulation of the latter |7j; and other recent remarkable progress both on it and on pending questions of 
the representability problem for 1-matrices [8||TT]. We return to the Ayers method in the next section. 



3 



here the use of a quasidensity or quasiprobability on phase space, which is nothing other 
than the reduced single-body Wigner function, whose relation to 1-electron space p(q) and 
momentum 7r(p) densities is immediate. 



Mathematically, Wigner functions 12 , 13 are equivalent to density matrices; and in this 
sense our approach, dubbed Wigner density functional theory (WDFT), is equivalent to 
1-density matrix functional theory. The beginnings of the latter go back to [14]; a surge of 
interest ensued, followed by a period of relative quiescence. More recently, this approach has 
seen a vigorous recent development centered around the notion of natural orbitals. However, 
WDFT can be worked out autonomously, calls for a different type of physical intuition, and 
readily lends itself to semi-classical treatments. Like one-body density matrix functional 
theory, it sits midway between ordinary DFT and the Coulson proposal of replacing the 
wave function by the 2-matrix fl5]. It stands to contribute to (and to gain from) both the 
Kohn and Coulson programs. 

Chief among our motivations is the success of the energetic program by Gill and coworkers 



see for instance 16, 17 — to attack the problem of correlation energies via the family of 



Wigner intracules. Their method is rooted in the fundamental observation by Rassolov 18 
that relative momentum is as important as relative position in determining interelectronic 
correlation; and this naturally calls for a phase space treatment. 

Since we assume from the reader only basic familiarity with Wigner quasiprobabilities, we 
do most things from scratch, starting from a proof of the Rayleigh-Ritz principle in phase 
space quantum mechanics. After this, we give the main theorems for energy functionals 
based on the single-body Wigner function, in parallel with standard treatments. There we 
are brief, since arguments of this type have become routine]^] 

The harmonium "atom" is treated exhaustively after that. Phase-space methods provide 
a fast lane towards an elegant and complete description of this system; it also serves as a 
training ground for employing the formalism. 

True atomic Wigner functions are considered next. Phase-space eigenfunctions which 
are not quasidensities are needed for that. We characterize them for Gaussian basis wave 
functions; this result seems to be new. We perform some numerical atomic computations 
with those eigenfunctions. 

Afterwards we (precisely define and) study natural Wigner orbitals. This is the heart of 



our subject. They come in handy to give a simpler proof than that of 19 of the "overbinding" 
property of the Miiller 1-density matrix functional for helium. Again turning to harmonium, 
we illustrate with them the workings of exact functionals. The Shull-Lowdin-Kutzelnigg 
series is analytically summed here, for the first time as far as we know. 

In an appendix, we place the Thomas-Fermi (TF) models within our framework. The 
purpose is mainly pedagogical, illustrating how they fit harmoniously within WDFT. An- 
other appendix elaborates on the characterization problem for Wigner functions. 

The fundamental ideas of the present endeavour are found already in |20] ; they appeared 



more developed in 21 . Wlodarz came upon the same idea by the mid-nineties 22 . But 
apparently his paper had an almost negligible impact. It is perhaps prudent to clarify that 
our approach is not quite in the same spirit as the (then recent) work on a phase space 



3 We consider fermionic systems only in this article. The parallel formalism for bosons will be considered 
elsewhere. 
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distribution corresponding to a given electronic density p, reported in the texts 23 24 . This 



had been conjured by Ghosh, Berkowitz and Parr 25,26 by means of heuristic "local ther- 



modynamics" arguments. The GBP Ansatz was put to good use 27 29 in the calculation 



of Compton profiles, exchange energies and corrections to the Thomas-Fermi-Dirac energy 
functional. Since the GBP distribution cannot be a true Wigner quasiprobability, its success 
should partially be attributed to the intrinsic strength of the phase space formalism. 

There might be an issue of the use of Wigner quasiprobabilities versus other correspon- 
dences between ordinary quantum mechanics and c-numbers on phase spacej^] Let it be said 
that, in the context of our quest, the mathematical advantages of Wigner functions are over- 
whelming. It would be difficult to develop natural orbital theory on phase space without the 



tracial property 31,32 underlying (uniquely) the Wigner-Moyal correspondence, that allows 
for an almost verbatim translation of the spectral resolution for quantum states. Also the 
metaplectic invariance of this correspondence apparently is essential for the reconstruction of 
the 2-body Wigner function from the 1-body Wigner function for harmonium in Section 
It lies behind generalized Wigner functions produced by fractional Fourier transforms 
as well. 

Last, but in no way least: this article seeks to continue and pay homage to outstanding 



7.6 



33 1, 



work of many years by Dahl and Springborg (see 34-36 and references given later) on 
atomic and molecular Wigner functions, by addressing several matters of principle. 



2 Classical Hamiltonians and Wigner functions 

The Wigner quasiprobability corresponding to a density matrix is given by 

iQn'iPij ■ ■ • > Pn'i ft> • • • > <>iv> ftj • • • > ^n) 
= ?r~ 3Af J D N (q x -Zi,sx,...,q N - z N , qv; q x + z u q[, . . . , q N + z N , q' N ) 
x e M?vzi+-+PN-ZN) d ^ d j N 

= n~ 3N J D N (1', . . . , N'; l", . . . , Ar") e i ^-(5i-5i)+-+^-(5^-5' JV )) 

x & (Qi ~ Wi + Qi)) "-S(q N - \{q N + q^)) dd 'i ■ ■ ■ d QN d Qi ■■■ d " 



N- 



(3) 



This is symmetric under particle exchange. The definition extends to transitions 



P^ n ^>' n (Qi7 ■■■■> Qn'iPh ■ ■ ■ iPn'i ft> • • • > Sn] ft> • • • , ftv) 



TT 



-3N 



*at(<Zi - zi, S!,...,q N - z N , q N )^/' N (q 1 + Zi, q[, . . . , q N + z N , q' N ) 



(4) 



We can regard Wigner quasiprobabilities as 2N x 2N matrices on spin space. When there 
is no risk of confusion the corresponding spinless quantity, obtained by tracing on the spin 



4 We still find most illuminating the treatment of the relations between phase space representations given 



by Cahill and Glauber 30 long ago 
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variables, is denoted by the same name. Then hermiticity of Dn in ^ translates into reality 
of Pn- The unit trace property translates into the normalization condition 

J PN(q 1 ,---,q N ;Pi,---,P N )dq 1 ...dq N dp 1 ...dp N = l. (5) 

Also, Wigner functions are square-summable and bounded continuous. By interpolation, 
Pn £ L P (M. 6N ) for all 2 < p < oo. Notice that the above integral at any point of phase space is 



the expected value of a Grossmann-Royer parity observable 37-39 ; those are Stratonovich 



Weyl (de)quantizers on Kirillov coadjoint orbits 31,32,40 . Since a parity operator cannot 



have expectation value greater than 1, we note that 

\P N (q 1 ,...,q N ;p 1 ,...,p N )\ < 7r~ 3N , (6) 

with negative values being possible. This bound cannot be reached in more than one point, 
and in view of (|5| it is often remarked that the support of a Wigner function is of bigger 
volume than (h/2) 3N . As a matter of fact, no square-summable Wigner function can have 
compact support in all of its variables, since then both the corresponding wavefunction and 
its Fourier transform would have compact support, which is impossible. 

The distinctive trait of the phase space formulation is that expected values are calculated 
as phase space averages. Hence the classical observables reenter the quantum-mechanical 
game. Let / be a real function of q x , . . . , q^lPi, ■ ■ ■ ,Pn ( a symbol) and let Q(f) be its 
quantized operator version (by the Weyl rule). We work with symbols symmetric in their 
arguments. Then we get: 



tr(Q(/)Djv) = / f{qi,---,qmPx,---,P N )PN{qi,---,qmPi,---,PN) dq 1 ...dq N dp 1 



>N- 



One must keep in mind that the symbol for the operator Dn is not P/v, but (2tt) Pn —'■ Wn- 
That is, Q(-P/v) = Dn/{27t) 3N . A Wigner quasiprobability behaves in almost every way like 
a probability density on phase space, except that in general is not nonnegative everywhere. 
It is not so easy, however, to recognize quantum state representatives among all functions 
on phase space: see Appendix [Bj 

2.1 Reduced Wigner functions 

Reduced Wigner quasiprobabilities, up to and including the single-particle quasidensity d 
on phase space, are obtained in the obvious way by integrating successively over groups of 
6 variables. In particular: 

4(?i,g 2 ;Pi>P2;^2;^,4) 

N(N-l) f .... . , ; , 

dq 3 . .. dq N dp 3 . . . dp N cfc 3 . . . ds N ; 
d{q\p;s\<?) = dx{q] := Jfzr\ J d 2 (q,q 2 ;p,p 2 ;q,q 2 ;q',q 2 ) dq 2 dp 2 dq 2 - (7) 
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By partial integration of d one gets the 1-electron density p and momentum density it: 



P(Q) = j ' d{q;p;q'S)dpdq; ?r(p) = J d(q;p;q;q) dqdq. (8) 

Also the form factors are easily expressed in terms of d\. Since dequantization and reduction 
commute, we may also use formulas analogous to (|3| for reduced density matrices. At each 
order the reduced Wigner distributions contain the same information as the corresponding 
reduced density matrices. It is convenient to have "inversion formulas" to recover the latter 
matrices from the Wigner functions. It transpires that 

r(l,2;l',2') = j ^^L^k^;^;^;^ 

7(1; 1') = I d(^;p;,;<;')eW'Up, (9) 



where T = t,N(N — 1)D 2 and 7 = ND\ as usually defined in quantum chemistry. (Note that 
d 2 corresponds to V rather than to D 2 and d\ to 7 rather than to Di.) 



One may usefully translate Coleman's theorems 



5] on representable 1-density matrices 



Let V]^ denote the set of 1-matrices 



into properties of d. This was done by Harriman 41 
representable by ([!]); it happens that V]^ CD 1 . In fact, d is the 1-body Wigner quasidensity 
corresponding to an element of T> X N if and only if 

0< J ' W(q;p)d(q;p)dqdp<2 (10) 

for any W which is the symbol of an arbitrary element of T> . Recall that d may be regarded 
as a 2 x 2 matrix in spin space. We can write for instance, 



and d(q;p) in (10) means the trace of d(q;p;q } <;') — or its scalar part, in a more cogent 
specification of d according to the behaviour of its components under rotation. Analogously 
d 2 may be regarded as a 4 x 4 matrix in spin space or as a direct sum of higher tensor 
representations of the rotation group. The Ayers trick — see [7| for instance — is easily 
reformulated here: the 1-body Wigner quasidensity belongs to T> l N if and only if 

h(q;p)d(q;p) dqdp > E [h; N] 

for every 1-body potential h(q;p), where E [h;N] denotes the ground state energy of iV 
fermions under such potential; likewise the 2-body Wigner quasidensity d 2 is iV-representable 
if and only if 

h2(q 1: q 2 ]Pi,P2)d2(q 1 ,q 2 ]Pi,P2)dq 1 dp 1 dq 2 dp 2 > E [h 2 ; N] 
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for any 2-body Hamiltonian ft, 2 - The "only if" part is clear, for otherwise the variational 
principle would be violated. The "if" part requires the Hahn-Banach theorem]^] 

The quasidensity corresponding to a Slater determinant (that is, an extremal element 



of V l N ) is the sum with equal weight 1 of the contributions of the occupied orbitals 35 

N 



i=i 

with the /j being (mutually orthogonal) pure-state Wigner functions. Note, as well, that in 
this case 

2N 

d(q;p) 2 dqdp = — -— . (11) 

In summary, quasidensities look like mixed-state Wigner quasiprobabilities. Typically, the 
latter are "less negative" than Wigner functions representing pure states. 

2.2 Spectral theorem and variational principle on phase space 



The spectral theorem for phase space quantum mechanics was given in 42 , 43] . Its formula- 
tion demands the concept of Moyal or twisted product x, indirectly defined by Q(/ x g) = 
Q(/)Q(flO- To alleviate notation in what follows we employ the convention q = (g 1; . . . , q N ), 
and similarly for the other variables. The twisted product is 



/ x 9(q;p) 



J /(</; P')g{q"; P") exp[2i(gp' - q'p + q'p" - q"p + q"p - qp")} 



Ir dq' dp' dq" dp" 



with the properties: 

f f * g(q;p)dqdp = J f(q;p)g(q;p)dqdp = J g x f(q;p)dqdp; fxg = gxj. (12) 

Theorem 1. Assume for simplicity a purely discrete nondegenerate energy spectrum E < 
Ei < ■■■ < E n < ■■■ . Then: 

• The solutions T nm , for n, m = 0, 1, 2, . . . , of the simultaneous eigenvalue equations: 

h x r nm E n r nm , Trim ^ h E m r nm , (13) 

form a doubly indexed orthogonal basis for the space of functions on phase space. These 
functions describe stationary states when n = m; when n ^ m they describe transitions 



between pairs of states. Note that T nm = T mn by (12). 

The sequence of those eigenvalues gives precisely the spectrum of H . 

With the normalization 

|2 dqdp _ 1 



nm 



(2< 



3N 



^Thus the axiom of choice. 
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we obtain 



r„ m x Tki = T n i 5 mk and the identity / T mn (q;p) f^J^ N = $mn, (14) 



consistently with (12) and (13). Clearly T nn / [2n) is the Wigner quasiprobability 



corresponding to the state of energy E n . 

• Also, the following closure relations hold: 

J2^nn(q;p) = 1; ^2^nm{q]p)T mn (q']p) = (2?r) 3iV 8(q - q 1 ) 5(p - p'). 

n m,n 

• Finally, we can write 

h ^ E n r nn . 

n 

Corollary 2. For any normalized state Pn and any Hamiltonian H , if Eq denotes the energy 
of the ground state corresponding to H , then: 

E < H(q;p)P N (q;p)dqdp. 



Equality is reached if and only if Pn = Pgs, the Wigner distribution for the ground state. 

The validity of this variational principle is of course not restricted to electronic systems. 

Proof. Any function F of the phase space variables may be expanded in a double series of 
eigentransitions corresponding to H . In particular, 

m,n 

where: 

P N (q;p)T nm (q;p) 



Thus 



(2tt) 3JV ' 

J H(q;p)P N (q;p) = J (H X P N )(q;p) *^ = J> 



where we have used (12), (14) and (15). Now take for F the symbol Wn of a pure state, so 
that Wn = Wn x Wn- An immediate calculation gives c nn = J2 m \ c mn\ 2 > 0. The rest is 
obvious. □ 
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3 The main theorems 

Let us invoke the classical Hamiltonian of an electronic system under the form: 

H — T + V ee + Vext, (16) 

where 

JV 1—12 1 N 

i=l i<j IVl i=l 

Here V ext denotes the external "potential" (e.g., due to the nucleus in an atom plus an 
external magnetic field). It is clear that the energy of a iV-electronic system is a functional 
of c?2- But we can prove more. 

Theorem 3. The many-body ground state of the system is determined by its 1-body quasi- 
density. 

Proof. Like in DFT, we argue by contradiction. Suppose that we have two different (by 
more than a constant) external potentials V eKt , V^ xt acting on our electronic system, with 
corresponding ground Wigner states P gs , P gs (assumed different) and respective ground state 
energies E , E' , whose associated 1-body quasidensity is the same. Then, by the variational 
principle on phase space of the previous section, we get 

E' = J H'P' gs dqdp< J H'P gs dqdp 

= J (H + VU - V cxt )P gs dqdp = E + J (\/' xt - V cxt )P gs dq dp. 

But 

J (Kxt - V cxt )P gs dqdp = J (V;' xt - V cxt )(q; p) d gs (q; p) dq dp. 
The same argument shows that 

E <E' + J (Kxt - V^)(q;p) d g MP) dqdp. 

and we obtain a contradiction. Thus d gs fixes P gs (and thereby the expected value of any 
observable in the ground state). □ 

Next we avoid l^-representability problems by a constrained-search definition of the en- 
ergy functional. 

Theorem 4. There exists a functional A of the electronic quasiprobability d, such that: 
E < J V(q;p)d(q;p)dqdp + J^d(q;p) dqdp + A[d] =: £[d\. 

Moreover, if d gs is the quasidensity corresponding to the ground state, then: 

f f \p\ 2 

E = V(q,p)d gs (q,p)dqdp+ / — d gs (q,p) dq dp + A(d gs ) = £(d gs ). 
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Proof. Let A[d] be 



min / P d (q;p)V ee (q;p)dqdp, 



(17) 



where Pd runs through every Wigner iV-distribution (representing mixed states in general) 
giving the fixed quasidensity d. Let P™ n be the one attaining this minimum. (Such P d form 
a compact convex set, with respect to a topology that makes Pj f PdV ee a continuous 
linear form, so its minimum will be attained.) The variational principle says that 

V(q;p)d(o;p)dqdp + J ^vr(p) dp + A[d] = J H(q;p)Pf n (q;p) dq dp > E . (18) 



E = H(q;p)P gs (q;p)dqdp< H(q;p)P^(q;p)dqdp 



In particular, 



This gives 

j P gs (VP)V cc (q;p)dqdp< J PZ n (q;p)V cc (q;p)dqdp. 
On the other hand, by definition, the reverse inequality holds: 

P dlT(q; P) v ee{q; P) dqdp< I P gs (q; p)V ee (q; p) dq dp. 



□ 



The minimization has been carried out in two stages. First we perform a search con- 



strained by the trial quasidensity d. In the second step, expression (18) is minimized. By 
Corollary 2 P gs = P™ n , which means that P gs can be obtained from P™ n directly even if 
the external potential is unknown. As in Levy's formulation of the Hohenberg-Kohn func- 
tional [44], A is universal, i.e., independent of V; nor is a ground state ^/-representability 
condition required. Thus our variational principle escapes the major problems of the original 
DFT one — see the discussion in |45| Chap. 33]. Systems with external potentials depending 
on momenta (as with orbital magnetism) fall into its purvey. Also, in the variation above we 
need not restrict to one-body Wigner quasiprobabilities corresponding to pure states; only 
finiteness conditions, such as that d should belong to the domain of the kinetic energy, 



T[d] 



d(q; p) dq dp < 00, 



(19) 



are implicit. Thus our apparent restriction to nondegenerate ground states is merely an 
inessential notational simplification: WDFT is an ensemble functional theory able to deal 
with degenerate ground states. For the same reason, A[d] is a convex functional, and there- 
fore S[d] too is convex. 

3.1 Exact requisites for the quasidensity functional 



In this subsection we collect a number of properties and conditions on A[d] within exact 
WDFT. 
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A Legendre transform variant of A exists: 

A L [d] =sup E[V ext ] - / (l\p\ 2 + V cxt (q;p))d(q;p) dqdp 



The supremum is taken over all possible choices of the external potential. The Legendre 
transform allows a reformulation of the previous variational theorems, in the spirit 



of 46 



Let V ext depend on a parameter A. Then we denote the Hamiltonian by H x . Consider 
the associated eigenvalue problem: 

H x xT x = E x T x = T x xH x , 

with T x being a normalized Wigner distribution. The assertion that 

dE x f dH x (q;p) 

-dX=] q\ q P (2°) 

is the Hellmann-Feynman theorem in phase space quantum mechanics. 
Proof. Indeed, E x = f H x (q;p)T x (q;p) dqdp implies 

dE x _ f dH x (q;p) ^ f^^^j^ , f u dT x {q;p) 
dX 

However, 



/ rA ^ 9 '^ d 1 d P + j H ^P) dTx ^ P ^ dqdp. 



/^(,; P ) S ^rf,* = /^( 9 ;p)x S ^ d , dp =£ J / Sr ^rf 5 <ip = o. 

The theorem follows. □ 
Now, let us write 

V cxt ,\(q;p) = V ext (q;p) + An • V^V ext (q; p) + Ar 2 • VpV ext (q;p) H 



for arbitrary vectors f\,r2- Thus, by (20) and homogeneity of space, 

d(q',p)(Vq + Vp)V cxt (q;p) dqdp = 0. 

Similarly, by isotropy of space, 

J d{q; p)(q x V 5 + p x Vp)V cxt (q; p) dqdp = 0. 

Only stationarity of the state is required. These results follow as well from the minimum 
principles of the previous section. 
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For stationary states, 

f \P\ 2 

E = - — d(q; p) dq dp 

= \\ I V( 3) d & P) dqdp+ [ d ll02^1l^l. dq 1 dq 2 dp x dp 2 

This is the virial theorem in phase space quantum mechanics; we consider here pure 
Coulomb systems. As a corollary, we get for the ground state: 

f \p\ 2 1 ( f 

E = - — d gs (q,p) dqdp = - f / V(q,p)d gs (q,p) dqdp + A(d gs ] 

We should not overlook that the minimization process takes place under the constraint 
J d(q;p)dqdp = N, which can be implemented by a Lagrange multiplier. Then one 
may minimize 

8[d] — /i — J d(q; p) dq dfij . 

The multiplier \i (with a minus sign) is an important physical parameter, called (Mul- 
liken's) electronegativity. Recall that for a neutral atom E (N — 1) — E (N) is the 
ionization potential, and Eq(N) — Eq(N + 1) is the electron affinity Their average 
constitutes a finite-difference approximation to the electronegativity^] 

Finally, we look at scaling. Matters are pretty satisfactory with d in this respect. Let 
A > be a scale factor. We scale the Wigner distribution by defining P\(q;p) '■— 
P\(Ag;A _1 p), another Wigner distribution, whose scaled quasidensity is d\(q;p) = 
d(Xq; A _1 p); yielding the scaled density p\(q) = A 3 p(Ag). One can show that d\ repre- 
sents a Wigner eigenstate of a Hamiltonian of the form T + XV ee + J2f=i ^^diq^Pi)- 
Now T[d\] is X 2 T[d]. As a consequence, 



A[d\] = XA[d}; and trivially 



dX 



= A[d}. (21) 

A=l 



The situation in the standard DFT approach with respect to scaling is much more 
involved. Denote Q[p] = T[p] + V ee [p], the universal Hohenberg-Kohn-Levy functional. 
The naive expectations T[p\] = X 2 T[p] and Ke[PA] = AKe[p] are both false: one can 
show that T[p x ] > X 2 T[p], V cc [p x ] < XV ee [p) for A < 1; whereas T[p x ] < X 2 T[p], 
V ee [px] > XVee[p) for A > 1. Nor is it possible to partition Q into two functionals in 
some other way with the desired behaviour |48| . 

We emphasize that these constraints on A[d] are valid for arbitrary quasidensities; 
this is why they are potentially useful. Explanation of the good behaviour in this 
regard of WDFT with respect to ordinary DFT lies obviously in the exactness of the 
kinetic energy functional in the former; whereas in the latter T[p] is a big unknown 
complicated functional, which "pollutes" the Coulomb energy. 



°The behaviour of these quantities is of current theoretical interest as a marker of the limitations of 



approximate functionals; see 47 and references therein. We cannot go into that here, however. 
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In summary, WDFT splits the problems of density functional theory into the (solved) 
characterization problem for 1-body Wigner functions and the determination problem for 
A[d], a functional both smaller in magnitude and less slippery in principle than that of 
Hohenberg-Kohn theory; exactness and simplicity of the kinetic energy functional commend 
our method. But one needs familiarity with the lore of Wigner quasiprobabilities. 



4 Getting used to Wigner quasiprobabilities 
4.1 Harmonium via the Wigner function 

In order to win intuition on the workings of A[d], it is good to study the WDFT functional 
in an analytically solvable problem. So we consider two fermions trapped in a harmonic 
potential well, which moreover couple to each other with a repulsive Hooke law force; this 
is the so-called harmonium, or Moshinsky atom 49 . The one-dimensional case has been 
treated on phase space in 50 . Introduce extracule and intracule coordinates, respectively 
given by 

#=(?i + <Z 2 )/2, 

with conjugate momenta 

P = p 1 +p 2} P = (Pi-p 2 )/ 2 - 
The classical Hamiltonian is given by 

p2 „2 

H = H R + H t = — + u 2 R 2 + p 2 + (co 2 - fc)-; 

the last term includes the electronic repulsion — Ax 2 /4 (we assume < k < u 2 ; obviously for 
k > uj 2 the repulsion between the particles is so strong that they cannot both remain in the 
well). The energy of the ground state is clearly given by E = \(u + ^u 2 -k). 

The corresponding Wigner function for the ground state factorizes into an extracule and 
an intracule phase space quasidensity: 

^gs(9i, fePuPri ft, ft! ft/ft") 

= ^(til 2 " Iit 2 )(ti42' " li't 2 ') ^ exp(-^) exp(- ^_ ). (22) 
Note the correct normalization 



J P SB (Qi, Q 2 ] V\3i\ ft, ^2; ft, ft) oil d2 



The electron interaction energy can be obtained from the intracule in (22): 



4vr 3 J r ^^ZF^k 



We have generalized the Wigner intracule, in the terminology of 17 , there valid only for 
non-interacting fermions in the harmonic well. In their notation, after integration over the 
angles it is given by 

w(u, v) = -u 2 v 2 e -^muy2 e -vV2v^^ (23) 

7T 
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Here v = 2p. By the way, the above type of calculation applies quite generally, not only for 
Hartree-Fock (Slater determinant) states, as declared in that paper. 

Now that we are at that, let us compute the relative-motion and centre-of-mass compo- 
nents of the kinetic energy Tq and the confinement energy. For the former we get, respectively, 



1 [ e -2H,,Vtt p 2^fi = MJ_Ji and JL [e- 2H «/"P 2 dPdR=^. 
7T d J 4 An 6 J I 



and the virial theorem is fulfilled, since En = 2Tn. For the latter 



^_ f e -2H</V^k r 2 ^ = and <J f e -2H R/ . R 2 ^ d p 3u.' 

4tt 3 J r A^/^P^k 7T 3 J 



In all, 



3a; 3a/oj 2 — k 3co 2 3k 3 /—= 

1 1 . = -leu + v or — k ), indeed. 

Since the value of the centre-of-mass energy is preordained, it should be clear that only (the 



two marginals of) expression (23) is employed. 



4.2 The Wigner 1-quasiprobability 



Also from (22), one gets 



PUvvfaPi,^) = ~e exp (~^~^i +p 2 + 2 Pi -A) - + <& + 2gi ■ g 2 )) 

x exp(- 2 ^__^ (^ + p 2 - 2 ft • p 2 ) - fa 2 + q 2 2 - 2q x • g 2 )J (24) 

Now we compute the 1-body phase space quasidensity for the ground state. One obtains: 



7T 



with marginal distributions 



(25) 



Pgs (q) = 2i- ^ ) e -2 gW - fe /(. + v^) ; (26) 
vr gs (p) = 2(un)- 3 / 2 ( % e -V/("+^). 



The normalization is now 



J d gM v\ ?; <0 dpds = J p g s(g) d v = J ^(p) 4p = 2 - 
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From the above we can easily recompute the kinetic and confinement energy parts: 



\p\ 2 7i gs (p) dp = 3 



UJ 



cxt,0 



\q\ 2 pg S (q)' 



UJ 2 + UO^/bJ 2 - 

Ayfuj 2 - k 



This yields Ke,o = —3k/A\/oj 2 — k, as we had obtained directly. 

Note that d gs (0; 0) < 27r~ 3 for k > 0; there is a telltale tassement of the Wigner quasiprob- 
ability with respect to what is typical for ground pure 1-particle states. Note also that 



d gs (q;p) 2 dqdp 



4 / 4wa/w 2 - k 



TP 



(u + V^r^ky 
Auj\/(jj 2 - k 



e 

3/2 



Aq 2 ujV^ 2 -k/(uj+Vuj 2 -k) -4p 2 / (ui+Vui 2 -k) 



< 



dq dp 



(27T) 



This is in agreement with (10). It is clear that this d gs cannot correspond to a Hartree-Fock 
(HF) state unless k — 0. 



We remark that the exact one-particle density matrix of 51 , Eq. (2-68)] and [52] for this 
problem is obtained from d gs simply by using the inversion formula (J9]). 

4.3 Pairs density 

Note that 



p(q x )p(q 2 ) = A(2uVi^k:/7i(uj + v^Tfc)) 3 e Z^S {9 * +q *\ 



On the other hand, by integrating out the momenta in (24) 

' lo^/iJ 2 — k^ " " 



P2@i,g 2 ) 

With spin components: 



7T- 



e 2 



\(w+\/w 2 -k)(q 2 +q 2 ) -{u)~yju) 2 -k) q r q 2 



(27) 



P2(l, 2) = -(ti4^ - I1T2XT1I2 " I1T2) P2(?i, q 2 ). 



Finally, we obtain 



without recourse to wavef unctions, density matrices, or the like. As was pointed out in 51 



besides the angular correlation in the pair distribution, favouring q 1 = —q 2 , which was to be 
expected, we see a contraction relative to the uncorrelated distribution. 



The above calculation can be organized in a better way, by integrating in (22) 

2H r 



7T U 



2H 



exp 

^V /2 p -2a; J R 2 



R 



cxp 



sfZF^k 



dPdp 



Vu 2 -k 
2^ 



3/2 



-^J2^fct 2 /2 



E(R)I(t), 



coincident with (27) 
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4.4 On the functional theory 



Now, according to the tenets of functional theories, V ext must be a functional of p gs or 7r gs ; 
depending on whether we employ ordinary DFT or Henderson's variant based on the mo- 

Dahl argues in 



mentum density 54 



50 



that (modulo an arbitrary constant Vo) one can 
indeed recover the potential, granted the harmonic form for it. But in general one would 
only be able to estimate the second derivative of the confining potential at the origin. So 
the Kohn-Hohenberg method remains nonconstructive, even for the harmonic interaction 
between the fermions. 

What does this example teach us about A[d] when V ee is of the harmonic form? One 
could be seduced by the following chain of reasoning. Note first that 

uj\ 3/2 { 2V^k 



Po:=Pgs(0) = 2( - 



+ Voo 2 -k 



UJ 



7T := 7Tgs(0) = 2 



UJ7T 



3/2 



2uj 



UJ + y/uj 2 — k 



3/2 



from which 



uj = uj(k; p /n ) 

is obtained by solving a simple algebraic equation. Therefore for such uj we have 



(28) 



T// -* uj(uj + \/uj 2 -k) po 



AVuj 2 - k 
and one could surmise that 



P 
~2 



UJ + \Juj 2 — k 7T 

A l0 S 



7T, 



^=^0- J V(q)p(q)dq- j 



-7T 



UJ 



UJ^ 



6 + / ir(p) Ioe 



7T 



7T 



dp + 



UJ 



^ r , p{q)\og — dq 

\/uj 2 - k J Po 



where in the last line all reference to the external potential has been banished. Even so, 
we have not obtained the universal Wigner functional for harmonic interparticle actions, 
because we have unduly restricted the variation defining A[d]. In spite of this, for confined 
two-particle systems the above formula doubtless constitutes a good approximation — in 



parallel to what was shown in 55 in the context of ordinary DFT. 

On a more narrow definition, restricting ourselves to harmonium, and given that we 
know the functional forms of the Wigner one-body quasidensity and the Wigner intracule, 
we certainly can determine the strength of the particle-particle interaction from d — it is 
enough to look at d gs (0; 0) — and thus the interaction energy. 



5 A Gaussian interlude 

Beyond the ease of calculations with Gaussians quasiprobabilities, there are pertinent rea- 



sons, from the use of Gaussian basis sets in standard quantum chemistry 56,57 , and from 
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entanglement theory 58 , that make it imperative to learn how to manipulate them. Also, 
phase space functions with negative regions can be reached by means of transitions between 
Gaussians. So we need to characterize those transitions. This is taken up in this section. 
We consider here s-type Gaussians centered at the origin; any others can be obtained by 



derivation — see the appendix in 35 — and translation. Not every Gaussian on phase space 



can represent a quantum state. Suppose that, with u = (q,p), we do have a Gaussian 

P F (u) = 7r" n (det F)a exp(-w*Fw), thus normalized by / P F (u) d 2n u = 1; (29) 
and suppose we want it to represent a pure state. Then F, beyond being positive definite, 



must be symplectic 59 . This will entail detF = 1. Recall that, if J denotes the canonical 
complex structure, 

r ( n 1, 

the matrix F is symplectic if FJF 1 = J. Any positive definite symplectic matrix can be 
factorized as F = S l S where S and S* are symplectic, too. In particular, such matrices F 
are symplectically congruent to the identity. The space of such Gaussians is of dimension 
n 2 + n. In our case, n = 3N. 

A Gaussian on phase space represents a mixed state if F is symplectically congruent to 

diag(Ai,...,A n ,Ai,...,A n ) 

with < Xi < 1. This was found in j60j. The space of such Gaussians is of dimension 
2n 2 + n. 

Given a symplectic (and positive definite) F, we can partition it into four n x n blocks: 

A B 
C D- 1 

Here the diagonal blocks are invertible, which is not always the case for general symplectic 
matrices; we choose the notation D~ x for convenience. Moreover A = A f , D = D t , C = E 1 '. 
Now, 

(b* D-^j J [b 1 D-^j = J im P lies BD = DB ^ AB * = BA > A = D + BDB 1 . 

We can know which wavefunction (also of Gaussian form) a Gaussian pure state comes from: 
Pp is the Wigner function corresponding to |^ r )(^ r | where, up to a constant phase factor, 

(detD) 1/4 

= —nW~ e M~k -Dq-y-BDq). 

In view of the above, D and BD are respectively positive and symmetric. For instance, the 
ground state wavefunction of the two-electron system considered in the previous subsection 
is given by 
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A similar formula holds in momentum space. 

It is of high interest to find the Wigner transition corresponding to |^ r i)(^ r 2| ) where 
and ^2 are Gaussians; they are bound to play an important role in calculations with 
Gaussian orbitals in quantum chemistry. All the information about \&2 is contained in 
their phase space partners P 1 = Pp 1 , P 2 = Pf 2 , so we can characterize the transitions from 
the parameters of the quadratic symplectic forms Fi, F 2 . 

Theorem 5. The Wigner transition P 12 (q;p) = P 12 (u) between two pure-state Gaussian 
quasidensities P 1 (m) and P 2 (u) given by (29) is of the form 



is a complex symmetric and symplectic matrix with positive definite real part, whose com- 
ponents are given by 



D 



.12 



B\2 
Al2 



(-§(£>! - D 2 ) + \{B x D l + B 2 D 2 ))D£, 

D 12 + B l2 D l2 B{ 2 . (31) 



Proof. The integral 



(* 2 I ^) = (detA) 1/4 (det^) 1 / 4 f ^i^+d^b^-b.d^ d 3N g 

converges absolutely, since the complex symmetric matrix P 12 = |(-Di+Z?2+*(Pi-Di — B 2 D 2 )) 
has positive definite real part \(D\ + D 2 ), and in particular Pi 2 is invertible. Note that 
B\ 2 Dyi = — — -D2 + i(B 1 Di + B 2 D 2 )) is also symmetric. 
The Wigner transition is P 12 = P*^, given by Q. Explicitly, 

p i2 (g;p) = ^i) 1/4 (detP 2 ) 1/4 



X / e -|(9-«)-(-Dl+tBiZ) 1 )(?-z) e -i(g+«)-(£)2-iB3£) 3 )(g+*) e 2i(p-a!) £N z 



_ (detPQV^detPa) 1 / 4 Dm f z . Dl2Zp 2iz.(p+B 12 n 12q ) «n 7 

~ jj-9N/2 J 

_ (detPQ 1 / 4 (detPa) 1 / 4 (p +BiaDl2? ). Dr i Wl2 ,) 

Here (detPi 2 ) 1//2 means the branch of the square root of det Pi 2 that is positive when Pi 2 
is positive definite. We have used the standard formula for a Gaussian integral, that is 
straightforward when P 12 is positive definite, and extends by analytic continuation to the 
case when the real part \{Di + D 2 ) is positive definite 61, Appendix A]. 

The formulas (31) show that the assembled matrix P 12 is indeed (complex) symplectic. 
The Wigner transition P 12 (u) is moreover square-summable, so that P12 itself has a positive 
definite real part. 

Whether or not the obvious reciprocal of this result holds is still an open question. □ 
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Notice in passing that when D\ = D 2 = D and Bi = B 2 = B, then Z} 12 = D also, and 
we of course recover 

P( q - p ) = n - 3N e ~<l-M e -p-D^ Pe -HB + Bt) P 

For the needs of quantum chemistry, it is largely enough to consider real wavefunctions of 
Gaussian form. If B X = B 2 = 0, then D 12 = \{D 1 + D 2 ), B 12 = - D 2 ){D X + D 2 ) _1 , 

and A 12 = D 1 (D 1 + D 2 y l D 2 + D 2 (D 1 + D 2 )~ 1 D 1 , so that 

P i2/ n \ _ 2 3JV/2 (detDi) 1 / 4 (det D 2 ) 1/A q . (Dl(Dl+D2 )-i D2+D2(Dl+D2 )-i Dl)q 
W>P) 1 ^(det(D l + D 2 )) 1 ft 

x e -2p-{D 1 +D 2 )- 1 p e 2iqi{D2-D 1 )(D 1 +D2)- 1 +(Di+D2)- 1 {D 2 -D 1 ))p 

The interesting new thing is the last factor: this will allow Wigner functions associated to 
the interference of P 1 and P 2 to become negative at some places. Also, depending on the 
overlap of P 1 with P 2 , they will exhibit damped oscillations when both p and q are large. A 
simple example with N = 1 will soon be useful. To 

| 9 9/ / cy-\ o \ 3/4 1 9 

P^ 2 {q-p) = —e- ai ' 2Q e~ p /ai ' 2 there corresponds V 12 (q) = ( — I e~^ 9 ; 

and then 

P 12 (g;p) = 7r- 3 (^) 3/ Ye-^ l((? - 2)2 e-5^^) 2 e 2 ^^ 



4ai«2 



3/4 



_3 / " tu: l ct 2 \ e -2(jr 2 «i«2/(ai+a2) <0 -2p 2 /(ai+a 2 ) <= 2j(a2~cn)/(oi+ a 2 )(?-p 

(ai + a 2 ) 2 

6 Atomic Wigner functions 
6.1 Gaussian approximations 

In Hartree units the Hamiltonian of a hydrogen-like atom is 

2 151 

The wave function for its ground state is well known: 

lM3) = -7=-e-*-; ^(p) = ^- , 72 _ Lr2 ,, 2 - (32) 
Therefore Pu(q;p) is given by 

-Z\q-z\ -Z\q+z\ -2ip-z _ f 1 1 r 2ig-g /oo\ 

^ 6 6 " 7T5 7 dZ |p-^|2 + Z 2 |^ +? |2 + Z 2 e " 

For a long time it had not been known how to compute these integrals in a closed analytical 



form, although in one dimension an analogous problem was solved 62 ; the geometrical 
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treatment via the Kustaanheimo-Stiefel transformation in 63 allows only to recover partial 



data for hydrogen- like atoms. A nearly closed analytical form is now given in 64 



However, one can approximate (32) by a sum of M normalized Gaussians, with real 
coefficients; thus 



M M 



Pg° @; P) = E cfP u (q; p) + J2 ^ (P ij + P ji ) p) ■ 

i=l i>j 

According to what we have seen in the previous section, it ensues that 

M 



c 2 e ~a i r 2 -p 2 /a l 



7T a *■ — ' 

i=l 

+ — > CicA- ^ttt e a * +a i e a > +a J cos — J —q-p . (34) 

71 V(«i + «i) / V + / 

It is clear that the (exact or approximate) result only depends on r,p,9, with 9 being the 
angle between q and p. It is then convenient to take r,p,9, together with three auxiliary 
angles, as the phase space variables. 

Let us briefly consider the case M = 1 first. This amounts to taking a trial state which 
is exact for an oscillator. For the energy: 



p 2 Z 



16 f /" 2 -ar 2 f P 



71 



The minimum is found at a opt = 16Z 2 /9tt, so the "equivalent oscillator" has frequency 
cu = 16Z 2 /9n precisely. It is equal to — 4Z 2 /3n, a pretty good shot at the correct —Z 2 /2, 
given the roughness of the approximation. At the origin pjy takes the maximum 
theoretical value 1/tt 3 . 

In order to visualize the quasiprobability, one considers the function F^\r;p) obtained 
by integrating over all angles and multiplying by r 2 p 2 . Its maximum for Z = 1 is found at 

(l/y/a^~t, y/a^) = (1.33,0.75) 

in Hartree units, that is 1.33 times the Bohr radius ao for the distance from the origin 
and 0.75 Oq 1 ^ for the momentum. Its contour map is given in 34, Figure 1]. This is a 
rather featureless everywhere positive function, that gives a poor idea of distribution of 
quasiprobabilities in the H-atom. 

Let us now try M = 2, allowing for oscillations in the _F ls function. In view of (34) we 

get 



Pif(q;p) = % g-^-f 2 /"! + £| e -« 2 r 2 - P 2 A* 2 



2 Cl c 2 f 4aia 2 \ 3/4 _2^ r 2 2 (2{a l -a 2 ) 

— — — e ai+Q 2 e "i+"2^ cos 

n 6 \{a 1 + a 2 ) 2 J V a i + a 2 
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There is no reason for c\ and c 2 to be of the same sign. However, it is intuitively clear 
that combinations of the same sign are energetically preferable; in particular, among the 
orthogonal combinations 



2 + 2((^) 3/4 



1 9 9/ 9 9/ 

tD -a 1 r z ~p z /a 1 _| e _Q ~ P ' Q2 



7T" 



7T" 



± 



4aia>2 
7r 3 \(ai + a 2 ) 5 



3/4 



2oiQ2 , 



/2( 



e ai+a 2 e "i+"2 cos 



\ a i + «2 



-9 -P 



P 9 will have the lower energy. We content ourselves with studying the radial phase space 
function 



Fg\r,p) = 8n 2 r 2 p 2 F P[f (r, p, 6) sin 9 d6 = —r 2 p 2 {c\ 

Jo 7T 



2 71 2( n 2^-a 1 r 2 -p 2 /a 1 + c 2 e -a 2 r 2 -p 2 /"2 



') 



32 Q Q 

_l r p C1C2 

7T 



4aiQ;2 \ 3/4 - 2011012 r 2 - — — » 2 . /2(ai - a 2 



(ai + a 2 ) s 



\ cti + a 2 



-pr , 



where jo(x) = (sinx)/x. We expect i-*i^(0; 0) = 7r 3 , implying the constraint 



1 = c x + c 2 + 2cic 2 



(«i + a 2 y 



3/4 



Using the formula 



/- P 2 ( 2(«1 — OL 2 ) \ 2 
e "i+«2 j mo dp 
V «1 + « 2 / 



cti + a 2 

one obtains the radial density of charge 
4 



27r( a 1 + a 2 f' 2 
16 



_ (ai~a 2 ) 2 2 
e 2(ai+a 2 ) r 



p(r) 



7T 



2 3/2 -air 2 , 2 3/2 -a 2 r 2 , o < \3/4 -^^rr 

qtt/ e 1 + c 2 a 2 e 2 + 2cic 2 (aia 2 ) 1 e 2 



«1+Q2 r 2 



/•oo 

with the reassuring normalization / p(r) dr = 1. 

Jo 

For the energy integrals, first there are the contributions 

'3a 2 



To these we add 



V n 12 



cti« 2 \ 3 ^ 4 3aiQ: 2 



4 V vr 



/ 



7T Z 



2tt 



3/2 



3^- 



(«ia 2 ) 7/4 



2(«i + a 2 ) \oti + a 2 7 (a 1 + a 2 ) 5/2 ' 

/aia 2 y /4 4ttZ 4Z(a!a 2 ) 3/4 
V 7T 2 / «i + «2 y/n(a 1 + a 2 y 
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reproducing the above result (35) when a\ = a 2 . 
Collecting our formulas, in the end we arrive at 



2/ 3a i ov Oil \ 2 (3a 2 nr7 ja 2 



E(a u a 2 ,c u c 2 ) = ci[-^-2Z\l-^ \+c 2 2 l^-2Z 



+ 2cic 2 ^v^- 



[a 



!« 2 ) 7 / 4 4Z( ai « 2 ) 3 / 4 



'(ai + a 2 ) 5 / 2 v^( a i + a %) / ' 
We recall that this expression is to be minimized with the constraint 



1 = c\ + c 2 + 2cic 2 



4a 1 a 2 x ' 



(«! + a 2 ) 2 / 

Analytically, it seems a hopeless task. Numerically, it is found that 

ci = 0.821230, c 2 = 0.274403; a x = 0.403059 Z 2 , a 2 = 2.664500 Z 2 . 

The corresponding energy is —0.485813 Z 2 au, a good estimate given the simplicity of our 
approach. 

We can also now minimize with the same method 

E(ax,a 2 ,a 3 ,c 1 ,c 2 ,c 3 ) 

2/ 3«i ov , 2 /3a; 2 /a 2 \ . 2 /^3a 3 ^ ja 3 




<{-t- 2Z ^) + ^-t- 2Z n) +ci \-f- 2Z ^ * 

V (ai + a 2 ) 5 / 2 V7r(ai + a 2 ) J \ (a x + a 3 ) 5 / 2 vM^i + a 3> J 



constrained by 



3/4 \3/4 , \3/4 



1 = C { + C 2 + C 3 + 2CiC 2 ? ■ + 2CiC 3 ■ ^ + 2c 2 c 3 ^ 

1 (ai + a 2 ) 2 y V.(ai + a 3 ) 2 y V( a 2 + a 3 ) 2 



Now one obtains 



ci = 0.647676, c 2 = 0.407884; c 3 = 0.070476; 



ai = 0.302753 Z 2 , a 2 = 1.362579 Z 2 , a 3 = 9.000725 Z 2 ; 
and E opt = -0.496979 Z 2 au. 

We judge this good accuracy for the ground states of H-like ions, showing the viability of 
the phase space approach; the rule of thumb "three Gaussian type orbitals for each Slater 



type orbital" 57 is fulfilled. Wigner transitions hold the key to serious computations with 
Gaussian basis sets in WDFT: they allow insight on the effects of "negative probability" 
regions for Wigner quasidensities at low computational cost. 
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Let us come back to the case M — 1, for the helium series. We choose for the quasidensity: 

2 „.„2 „2/„, ., . ,_ s 2 /" _^2_„2/„, 3 _ 2a 3 / 2 _ ar 2 



d(q;p) = 2 e - Q r2- p 2/« 5 SQ that . = 2 /" e -c^-p>/a^ 

For the energy, first of all we get for the noninteracting part 

Bnl(a) = __4zJ-, 



just by multiplying the result in (35) by 2. The minimum of just this for the neutral ion 
Z = 2 would be found at a = 64/97T and is equal to — 32/37T. With the assumption of a 
singlet state there is no exchange energy, and the Coulomb electronic interaction integral is 
easily taken care of: 

4 J \q~tf\ A J t V 7T 

For Z = 2 this is 4a/2 times smaller in absolute value than the nucleus-electron energy, a 
smallish^ but roughly satisfactory ratio. We thus get 

E tot = — -(AZ-V2)J-; 

the minimum is now found at a = (16Z 2 + 2 — 8Z\/2)/9ir and is equal to 

4ZV2 -1-8Z 2 
3^ ' 

Thus for helium we get a ~ 1.53 and E ~ — 2.3au: far above the true energy, although 
hardly worse than the comparable result for hydrogen. Also, we already know from (11) the 
necessary equality 

d{q\p\q,qfdqdpdq = J d(q; p) 2 da dp = — J e ~ 2ar2 - 2 P 2 / a dq dp = . 

What we have just done coincides with the HF calculation and discussion of the corre- 
sponding Wigner intracules in 17, Sect. 7.1], in which the trial state is exact for a pair of 
uncoupled oscillators. In this reference the Wigner intracules for noninteracting fermions 
in a harmonic well is computed in closed form until N = 8 and they are ". . . surprisingly 
similar to those of qualitatively analogous atoms". Now, in such a context it might seem 
tempting to recruit to the cause the exact ground state for interacting harmonium, given 
that the intracule formula for the interelectronic Coulomb energy is very simple: 



Then E tot (a, k) = E n i(a, k) + E^ e (a, k), where the first part of the energy has the form 



, 4aVa 2 - k \ 3/2 / 3(a + \/ a 2 - k) „ / 2ayfa 2 - k 
E ni (a,k)=[ / ■ _ ^ >--AZJ 



[a + y/a 2 - k) 2 J \ 4 y 7r(a + V« 2 - k) 

However, numerical calculation show only a marginal improvement in the energy. 



7 This ratio is equal to 6.4 for the Kcllner model of He, and is bound to be larger for the "true" model. 
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Figure 1: Contour curves in hartrees for the Wigner function for the beryllium atom 
6.2 Real atoms 



For the H-atom, one may consult now 34, Figure 5] for F ls (r,p), drawn from a very good 
set of Gaussians with M = 10. The function there attains its maximum value at (1.30, 0.6). 
There is an infinite region of damped oscillations going into negative value regions, starting 
with a nodal curve going through r = 0.5 ao at p approximately equal to 4; through r = 1 
at p ~ 2.3; through r = 2 at p ~ 1.4; through r = 3 at p ~ 1; through r = 4 at p ~ 0.8. The 
amplitudes are small (_F ls (1.8, 1.8) = —0.0047); but oscillations are definitively there. 



Images for closed-shell atoms, based on Hartree-Fock configurations, are given in 35 
Experience with atomic ground-state Wigner functions allows one to reach the conclusion 
that the phase space region supporting an atomic Wigner quasidensity separates typically 
into three rough subregions. In the inner region the function Fi s mostly takes positive 
values; but it may take negative values, due to complicated interplays among orbitals: the 
one-body Wigner distribution for the ground state of neon exhibited in [35] is a case in point, 
because of the weight of the 2p orbital. One sees that in the dominant middle region, this 
function takes large positive values; but negative values may also appear, due to entanglement 
between electron pairs. In both these regions the distribution barely oscillates. In the edge 
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or "Airy" region we find an oscillatory decay regime. This decay of the Wigner distribution 
has been rigorously proved to be generic for systems with exponentially decaying (with linear 



exponents) states 62 , such as those of atoms. As it turns out, the middle region corresponds 
to the region of the atom where the semiclassical approximation is reasonable; and the nodal 
curve by the frontier between the middle and the outer regions reproduces surprisingly well 
the border of the occupied region in the TF model of the atom — even for the H-atom! 

In the inner region the Wigner function definitely departs from the Thomas-Fermi 
Ansatz, in that it is always bounded; also, as soon as N > 2, one finds "holes" of nega- 
tive quasiprobability not far from the origin in phase space. In Figure [I] (taken from |35|) 
the hole produced by the interference of the Is and 2s pairs is clearly visible. 

On the strength of all this work we realize as well that WDFT provides an easily visual- 
izable bridge between the TF and quantum approaches. We formalize this in Appendix [A|p| 



7 Natural Wigner orbitals 
7.1 Preliminaries 

Now we gear up for a new approach to A[d]. Generally speaking, this is equivalent to looking 
for a functional d2[d], or for a functional P2[d] (one does not need g?2 to compute the Coulomb 
repulsion, its diagonal is enough); or even a functional I[d], with / the position intracule 
(for which we presume there is no representability problem), would be enough. In view 



of the Schmidt theorems on best approximation 51,66 , it is natural to use the spectral 
representation of d(q;p;q;q'). Coleman's representability theorems [5] can be construed as 
implying the natural expansion 



d(q; p; q, = ^ «i fi(q; p; q') 



i=i 



where the occupation numbers fulfil < n» < 1 and Y2 n i = N\ we order them by 
descending size. A state with rij = 0, 1 is a pinned, extremal or HF state; we already 



discussed them in subsection 2.1 Examples are known of interacting systems for which 
0, 1 for some of the natural orbitals occurj^] For Coulomb systems typically the first rij for 
1 < i < N are close to 1, corresponding to a state close to the best HF state; and the others 
are small. A proof of infinitude of non-zero occupation numbers in this context has been 



claimed in 67 . There are tantalizing cases, however, where the best numerical computations 



stubbornly yield reduced states of finite rank — see 68 for the first excited state of beryllium. 

In the above sum, i carries both spatial and spin indices. To have pure spin eigenstates, 
the non-diagonal spin blocks must be zero; moreover, in a spin-compensated, closed-shell 
situation the diagonal blocks are equal. Then N must be even. In this case the spinless 
quasiprobability is of the form 



8 The nodal structure of the electronic Wigner function for molecules has been investigated in 
9 This point was clarified to us by the referee. 
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where the are still normalized to 1, but now < z/j < 2. Also, the fi verify 

d x fi = fi x d = Uifi. 



(36) 



We call these fi the natural Wigner orbitals (NWO) — natural Wigner spatial orbitals would 
be more precise. We shall also need the natural Wigner spatial transitions, denoted fij. We 
know that in principle they can be found from the fi and, for d spinless, they satisfy 



d x fy = Uif 



fij x d — Vjfij- 



7.2 The Miiller functional 

For d a HF state, the general form of A[d] for an electronic interaction / is well known: 



A[d] 



p(<ZiM92 )/(?!> ?a) dqt dq 2 - / (7(1, 2)| 2 f(q ly q 2 ) dl dl 



D(p,p)-X( 7 ), 



(37) 



where p and 7 are known functionals of d; the term X{^) is the exchange functional. If one 



uses (37) for arbitrary allowed quasidensities, one gets a functional for the interaction energy, 
proceeding through the pair density p<i[d], denoted *4. HF . Adding the kinetic and external 
energy functionals, we get a functional for the total energy, denoted £ HF [<i]. Provided / is 
positive semidefinite, its minimum is always reached in a HF state [53]. This means that 
£ HF [<i] renders an upper bound for the energy of the system; and that Slater determinants 
literally do no more than scratch the energy surface. The properties of £ HF [<i] are not very 
good, besides. Clearly it is not convex. Also, it does not respect a basic sum rule in general 
(this will be recalled in the next subsection). The difference between the minimum for the 
total energy attained by use of this expression and the "true" binding energy is generally 
called the correlation energy] it is obviously always negative. 

The problem of finding p2[d] has been considered in the context of one-body density- 
matrix functional theory. In a remarkable paper Miiller in 1984 [69] proposed an approxi- 



mate formula for D 2 amounting, in our context, in the notation of (37), to the alternative 
functional: 

A M [d} = D(p,p)-X(^). (38) 

Note that because of the equivalence of operator and Moyal product algebras, ^7 makes 
sense (meaning that ^7 x ^7 = 7), provided one treads carefully. The concept is akin to 



Wlodarz's "phase space wave function" 22,70 



After a period of some obscurity, the Miiller functional was rediscovered 71 , 72 and 



seems to be still in fashion 19 . The Miiller functional is indeed convex 19 . Clearly 



*4 M [<i] = *4. HF [<i] on the subset of extremal quasidensities. However, in general ^4 M [<i] < 
*4 HF [<i], and *4. M [<i] actually tends to give lower bounds than the true values. To prove this 
in general would be important]^] It is only known for sure for N = 1,2 as yet [19] ; and, 
as perhaps could be expected, the proof in this last cited paper is quite complicated. We 
show within this section a much simplified proof by use of the exact functional, formulated 
in terms of NWOs. 
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7.3 General NWO functional theory 

Let the spectral representation of d be given, and denote 



K = J + v {<ij) dqdpdqdq. 



The integrals are assumed to be finite. Then the energy is given by J2i n ihi + Kc^]- 
We want to express the electron-electron repulsion ^[^2] in the NWO representation. The 
expansion of d 2 in any orthonormal basis of eigentransitions, in particular the natural Wigner 
eigentransition basis, is denoted 

d a (l,2) = ^D2/«(l)/y(2) ) 

ijkl 

including spin indices]^] Note the symmetries 

DZ = D%*, D% = -D^ D% = -Dl D% = Djt; (39) 
so one can rewrite the above expansion as 

E DZ{f ^(1)^(2) + &(l)/*(2) - Ml)fu(2) - /«(l)/ w (2)). 

i<j, fc<i 

Notice that in this language, the HF functional corresponds to taking, in the natural basis, 

In principle d 2 has sixteen spin blocks, but, as a consequence of requiring pure spin states, 
only six differ from zero: 

A TT A 44- A 44 A 44 A 44 44 
^t-l-' 24-f 2i-f-) "2-f-i) 

and only three of those are independent. With this notation, Ke[^2] is given by 

YM 1 pg) K;? + <f| + <;? + 

where 



I9i -92 1 

with purely spatial indices. Note that the other two nonzero blocks cannot contribute to 
this Coulomb integral. At this point, it is convenient to introduce spinless density matrices 
by tracing: 

v v - ript + TXpj., .- i^ rt)tt + 1^4,4 + L> P ,Mt + D rm , (4UJ 

and write 

E = Y^h P + Y D r^ rt \P^- 

p pqrt 



11 Mutatis mutandis, we borrow the notation in the excellent reference 
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here. 
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Note now that the sums: 




ijkl ' ijkl ij 



must be N(N — l)/2, the number of electron pairs, while 



N 2 



En 2 



JV(JV-l) 
2 



> 



2 



2 



ijkl 



unless all n, are or 1, that is, the pure HF case. This is the sum rule we alluded to in the 
previous subsection. Let us baptize our ignorance 



this cumulant deserves perhaps to be called the correlation matrix. We now observe that, 
while certainly J p(l)p(2) all al2 = 0(N 2 ), we have 



This kind of observation is useful in Thomas-Fermi theory. 
7.4 NWO for ground states of two-electron systems 

For two-electron atoms we can do much better. Let us invoke invariance of the Hamiltonian 
under time inversion — the latter is represented by an antiunitary operator related to spinor 
conjugation that we need not write out. This assumption is not essential, but simplifies 
matters, and moreover holds in most cases of interest; it entails that the eigenstate wave 
functions are real. It is most instructive to start from a general basis of eigentransitions and 
construct the natural basis out of it. This is equivalent to recasting the results of the classic 
work [76] in our language. 

So let an orthonormal basis {f nm } for single-body functions on phase space be given, 
arbitrary except for the properties of Theorem [l] Consider singlet states. We want to 
expand a normalized singlet 2-Wigner function P 2 = d% in terms of the f nm . Its spatial part 
must be symmetric under exchange of coordinates; thus P 2 is of the form 





(41) 



^(gi^Pi^Shft^sa) = 2(^1^2 — iit 2 )(ti'i 2 ' - li>tv)f(qi,fePiiP2)i ( 42 ) 



where 



f(3 1 ,q 2 ;p 1 ,p 2 ) = - ^^^[/^(^ftj/ri^jft) + fri(qi,Pi)fmk(q 2 ;P 2 ) 



klmr 



+ fml(q 1 ;P 1 )frk(q 2 ;P2) + /r*(?i;Pl)/ml(?2;^2)]> ( 43 ) 
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with Cki = Cik real and J2h = 1- Each group of terms is real in view of f m k = fkm and 
the normalization ensures that 

J P 2(Qi, ^2\V\3i\ ft, ft, ft, ft) dq t dq 2 dp 1 dp 2 dq 2 = 1 



-M* I -4-1 I -t* -i 

in (42). It is clear that in the singlet case c^jj = c?2^ = and that c^jj = c^j = 2/ ■ 



The corresponding 1-Wigner distribution is obtained at once by integration, 

d(g- p- ? , o = (tf + II') £ 7^/^(9; p) = ( Ekm lm i mM ® „ °, „ ^ 

V Lfc m 7mfcimfc(<?;p), 



where 7 := [y m k\ = [X^z CmiCa] is a symmetric (and positive semidefmite) matrix. The 
matrix C determines a real quadratic form that can be made diagonal by an orthogonal 
change of basis O = [Oi m ], 



C = OcO\ implying 7 = O c 2 0\ 

where c = [c r S rs \. These elements c r are real, but note that some of them may be negative. 
Let us now make the definition: 

~Xrp • ^ OmrfmkOkpi SO that fkm ^ O kpX.'prO mr • 
mk rp 

The diagonal Xrr = Xr are real: %. = J2 m k O m rfkmO k r = Xr, and more generally Xr~s = Xsr] 
and it is easy to check that 

Xi X \j ^ ^ OmifmkOki X ^ ^ Olj fl s O s j ^ O m if ms O s j Sij Xi ^iji 
mk Is ms 

and that, with denoting either of the two nonvanishing spin components of d, 

d* ^ ^ ^Ymkfmk ^ ^ O mr '~y rr ikO kp Xpr ^ Cj, 5rpXp)* ■ ^ ^ ^-rX?" 
fcm rp km rp r 

(on writing n r := c 2 ), so the x% solve the simultaneous equations 

d* x Xi = Xi x = 



which is (36) up to a factor. Furthermore, J^nj = tr(c 2 ) = tr(C 2 ) = 1. The rij are the 
"occupation numbers" for "natural Wigner orbitals" Xu we have recovered the Coleman 
results in this case. It remains to compute 



^2 C mrC k l fmk(qi,Pl)frl($2>P2) = c vC t O mv O rv O kt Oi t fmktfliPl) J riGfa^a) 

vt,klmr 

^2c v c t Xvt(qi,Pi)xvt(q 2 ;P2)> 



kl.mr vt,klmr 



vt 
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and similarly for the other terms in (43). 

Since c r = ±\/7V, in order to obtain d 2 [d] and thus A[d] one needs "only" to determine 

Recall that in Q, in view of (42), the terms Z^fjf and D^jf 



an infinite number of signs 



12 



vanish. We see that D v r \ 



c r c p 5 r t S P9 in the new basis. We use the notations 



Xpr ( 



x P r(q,p) 



(rr | pp) 



Xpr(qi)X P r(q 2 



2- 



(44) 



borrowing now the eigentransitions \ vr in the expansion of d%. Note that these terms arise 
from correlation between electrons with opposite spin. They are real: in view of Q, the 
integrals over the momenta yield real functions since the corresponding wave functions are 

and L rp = L pr , too. In the weak correlation regime, when 



real. In particular Xpr(q) = Xr P { 
there is a dominant state close to the best HF state j^jLowdin and Shull [76 77 empirically 
found long ago that, if conventionally c\ was taken equal to +y/rii, then (most of) the other 
signs were negative. This gives rise to: 



g?2l(1,2) = (spin factor) x 



i(9i)xi(?2) - Yl V n i n p(xpi(qi)xpi(q2 



p>2 



A L [d] :-- 



+ Y y/ n r n pXrp(qi)Xrp(q 2 ) 
p,r>2 

n\L n - 2 22 \/ n i n P l ip + V n ^ n p L r P - 

p>2 p,r>2 



xi P {qi)xip{q 2 )) 

(45) 
(46) 



In this connection the work of Kutzelnigg 78 was decisive as well. As mentioned before, the 



natural orbital construction guarantees term- by-term the most rapid approximation to d 2 - 
The case where the sum stops merely at p = 2 accounts with remarkable accuracy for a 
good fraction of the radial correlation energy [76]. We might call these the Shull-Lowdin- 
Kutzelnigg functionals (SLK functionals, for short). For the comparison with the Miiller 
functional in the next subsection, the configuration of signs turns out to be irrelevant, and 



generally we call the analogue of (46) with any sign choice a SLK functional. 

The case for the (strict) SLK functional has been argued for mathematically as follows. 
By the Rayleigh principle, the quadratic form representing the energy has a minimum eigen- 
value when (ci c 2 • • • )* is an eigenvector, that is, 









>2nihi + Ln L 12 


■A 




( Cl \ 




c 2 




Li 2 2n 2 h 2 + L 22 






c 2 




w 




V : : 


■J 




W 



Therefore, for all r: 



Eq = 2n r h r + L rr + ^ L 
p¥=r 



rp 



12 This amounts to an "inversion" of the Schmidt decomposition 66 popular in entanglement theory. 
13 In which precise sense close, and how close, was discussed in 75 . 
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Only the two-body part of the Hamiltonian contributes to the off-diagonal terms. Well- 
known properties of the Coulomb potential ensure that these are all positive. Indeed, the 
Coulomb potential is positive definite in that its Fourier transform is a positive distribution. 
Or one can use its integral representation [791 Chap. 5]: 

7T 3 r i i 



dfz. 

I9l-92l J \ Z ~ <ll\\ Z ~ <l2\ 

Thus, whenever n\ (that is to say c\) is dominant, to minimize Eq we must put C2 negative, 
and probably C3, C4 . . . as well. We remark that in the weak correlation regime the energy 
matrix is diagonally dominant, and diagonalization of the quadratic form proceeds without 
obstacle; all pivots will be negative. Diagonal dominance ensures that the energy matrix is 



negative definite 80 as well. Nevertheless, while this reasoning shows that there must be 



some negative signs, it does not prove that all signs beyond that of c\ are negative. 



A different argument to the same effect was put forward in 81 . We summarize it next. 
Consider the gradient of the energy, 

4/i^>c^> I 2 ^^^^ CpL/ipp 2c^.A, 



da 



r 



where A is a Lagrange multiplier, at the "Hartree-Fock point", defined by c\ = 1, C2 = 
C3 = • • • = 0. For r > 1 only the second term on the right hand side contributes. Thus 
the only way to obtain a negative gradient is for c r to become negative in the minimization 



process. Also in 81 it is remarked that for the negative hydrogen ion H~, wherein the 
largest occupation number differs significantly from 1, the sign rule numerically holds true. 
Again, the argument is persuasive, but not conclusive, since a moment's reflection shows 
that it refers to the approximation from HF states rather than to the exact state. 

We now list some properties and traits of the SLK functional for singlet states. 

1. The sum rule J p2h(l, 2) dl d2 := f c?2l(<7i> Qz'-iPiiPz) ^<Zi dpi d% = 1 is (of course) 
fulfilled. 

2. The pair coincidence probability is a perfect square for any rank (and actually any 
choice of signs). 

3. For pinned states, it reproduces the results of the HF functional. 

4. The correlation energy is negative. 

5. The SLK functional scales linearly: 

A[d x ) = XA[d], 



with d\ as in subsection 3^ So the correlation energy functional must scale linearly, 
too. 



6. The SLK functional satisfies known constraints from the -D2-representability theory [5j. 

The previous one is a typical closed-shell configuration, and for such both spin compo- 
nents are completely identical. However, the fact of the matter is that occupation numbers 
always are evenly degenerate for any two-electron state jH]. 
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7.5 Shull— Lowdin— Kutzelnigg versus Miiller 



Let us write the Miiller functional (38) in terms of NW orbitals: 

P2u(qx,q 2 ) ■= \ Yn j n k x j {q l )Xk{q 2 ) ~ 9 $3 V n J n * Xjk{qi)Xkj{q 2 ) ■ 



2^ J ^^^"^ 2 

jk jk 



The integral over all coordinates comes out right as 1 (see below). But clearly d 2 M does 
not possess the antisymmetry properties required in ((391). As recognized in the original ar- 



ticle 69 , the corresponding pair density /?2m(<?, q) can take negative values for some states. 
Thus it is not surprising that the Miiller functional tends to give lower bounds than the true 
values; when applied formally to the hydrogen atom, this is the case. Unphysical probabil- 
ities together with the overbinding property clearly indicate that it does not correspond to 
physically realizable states. 

Reference 19 amounts to a long analysis of the Miiller functional, crowned by a proof 
that it gives lower energy values than the true values for helium. Now we deliver the promised 
simpler (as well as more informative) proof of that theorem. 

Theorem 6. For the isoelectronic helium-like series, the Miiller functional A^ld] is a lower 
bound to quantum mechanics. 



Proof. We simply compare directly the Miiller and SLK functionals. In the notation of (44), 
with Pi(q) = j Xiiq-iP) dp too, the pair density coming from the SLK functional is 

P2l(mi, Q2) = Y n ipi(qi)pi(q 2 ) ~ 2 Y1 V n ^ n j xij(qi)xij(q 2 ) + Y Xij(qi)Xij(q 2 ), 

normalized by / p 2 h(qi, q 2 ) dq x dq 2 = 2_^ ni = 1- (47) 
J i 

For the same object in our (two-electron, closed-shell) context, the Miiller functional gives 

p2M(qi,q 2 ) = ^p(qi)p(q 2 ) - J^p^iM^) - ^VwhXij (qi)Xij(q 2 ) ■ (48) 

i ij!=j 

It is perhaps not obvious which of the contributions \p(q\)p(q 2 ) or 2 £\ niPi(q^)pi(q 2 ) 
will be larger. However, on using nf = rii — J2j^i n i n jy we compute the "defect": 

2^^Pi(gi)pi(g 2 ) - \p(qi)p(q 2 ) = 2 Y niPi ^ Pi ^ ~ 2 {^2 n ^ d S) (2 n ^"@ 2 

i i i j 

This must still be integrated with 1/1^ — q 2 \; then we see that the positivity of the Coulomb 
potential alluded to in the previous subsection does the job, establishing that *4 M [<i] < 
A L [d]. □ 
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Remark 1 For a chemists' chemist, the proof is surely over. A mathematically minded one, 
however, would argue that we still ought to provide for convergence of the series representing 
the difference in energy. To guarantee that, consider the expression 



5 = E 



riiUj 



\Qi ~ 92 1 



dq 1 dq 



2- 



(49) 



where Fij := pi — pj. The Hardy-Littlewood-Sobolev inequality |82j Sect. 2.3] gives, for a 
suitable (sharp) constant C*i, the bound 



S < C\ riiiij \\Fi 



'ill 6/5 — Cl / ] n i n 3 \\Pi~ Pi 1 1 6/5 



< 2C 1 ^TwSPi\\*/5+\\PA\l/fd= 4C l'52 n i\\Pi 



|2 

1 6/5 " 



i>l 



As a marginal of a Wigner function, each pi is positive and lies in L X (]R 3 ), with = 1. 
However, we need to confirm that each pi(q) belongs to L 6//5 (IR 3 ). Finiteness of the kinetic 



energy produces the estimate we need for that. As suggested already by Moyal 13 , from 
formula (|3j) for (spinless) momentum states we obtain: 



P 

^ > / -jfi&P) 



2vr 3 



p 2 'j(p;p")e 1 ^'^ p ^ 8(p — \{p + p")) dqdpdp dp' 



-\ J (V 5 -V^) 2 7 (f;x')| 2=2 , =? dg. 



On account of the natural orbital expansion j(x;x) = ^ni^a;)*^^'), with pi = the 
right hand side can be written as the integral over q of 



i>l 



ri; 



V 2 p, 



and its finiteness entails 



y^^l|V(A/p7)||2 = j iy^/pi) 2 dq < oo. 

i i>l J 



The Sobolev inequality (for instance, see 83, Sect. 10.3]) can now be invoked to show 
that ^j^iPi(g) lies in L 3 (IR 3 ): 



for a suitable constant Ci- In particular, each pi lies in both L X (IR 3 ) and L 3 (IR 3 ), and by 
interpolation, it also lies in L 6 / 5 (IR 3 ). The related Holder inequality is 



Pi II 6/5 < HPilll llPilb' > 



3/4 1 1 1 1 1/4 
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and this allows us to complete the estimate of (|49j): 

i>\ %>x 
< 2C 1 ^n i (||p i ||? + ||pi|| 3 ) =2C l (l + ^2n i \\ Pi \\\ 

i>l ^ i>l ' 

since we have already shown that the last sum is finite. 

Remark 2 In the literature there is at least one attempt |84] to compare both functionals, 



predating 19 , by the way. According to 84 , two apparent differences should account for 



the Miiller functional's tendency to overcorrelate. 

1. The Miiller functional has more negative signs. Indeed, for the SLK functional(s) it 
remains that some weakly occupied orbitals contribute + signs. 

2. As we have seen in our proof of the theorem, the correct expression ^\ TijPjtf^Pjfy) 
is replaced by the approximate ^p(q 1 )p(q 2 ) — Y^j n jPj{di)pj{q 2 ) — the latter approxi- 
mation being regarded as largely valid in (84). Note that this partition explains why 
the sum rule is preserved by the Miiller functional: on integration, the crossed terms 
give zero contribution, and 

2 Yl n i J pMi)pM2) dd i dd 2 = 2 = \ J ptfMtiz) d Qi dq 2 - 

Curiously, the extra minus signs in Miiller's functional do not seem to play a role. 
7.6 Harmonium test for the Shull— Lowdin— Kutzelnigg functional 



To the best of our knowledge, the SLK series (45) and (46) have never been analytically 
summed. Here through the quantum phase space formalism we show that the harmonium 
model provides a first example of such summation, on the way deciding the sign dilemma 
for it. 

The 2-representability problem is the same as for real atoms, since it only involves the 
kinematics of fermions. On the other hand, the positive definiteness property of the Coulomb 
potential is lost; this and the confining nature of the one-body potential make for a peculiar 
determination of signs. We ought to effect a two-step procedure: 

• To expand the 1-quasiprobability d in natural Wigner orbitals. 

• Then to sum the SLK series to see (whether and) how the known expressions for d 2 
and the energy are recovered. 



The difficulties are just of a technical nature. We had the 1-quasidensity (|25|): 

2 (Aupf/ 2 

7T 3 (U + pY 



d gs {q,p- q, = (spin term) x 4r ^^[l e ~ 2 ^ q2+p2)/ ^\ (50) 
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where now we write /i := a/co> 2 — k. Since the three variables separate cleanly, we can drop 
the spin term, omit the normalization factor, and go to the one-dimensional case 

dJu) := dJq,p) := !V§| e - 2 ( UH ^V)/Kri 

7T OJ + /i 

Returning to the three-dimensional d gs is then just a matter of notation. 

From (22), for the 2-particle density g?2*, with the spin term omitted and reduced from 
three dimensions to one, one gleans: 

d 2 *( Ul ,u 2 ) = L e -™*/» e -2Br/r (51) 
= ^ ex p(-|(9i + <? 2 ) 2 - -g 2 ) 2 - ^-(pi +P2) 2 - ^-(pi -p 2 ) 2 ^- 

Notice in passing that this is indeed a pure state, since here e~ u ' Fu = e~ q ' Aq ~ p ' A lp , where 

+ oj — fi\ _ 1 /or 1 + UT 1 - /i -1 

~2W-/i uj + jj,)' ' ~ 2 la; -1 - w" 1 + /x -1 



Formula (51) is what we need to recover. Since mathematically (50) is a mixed state (in fact 
a maximally mixed one, as we shall see, for the given value of the parameters), the problem 
is not unlike mending a broken egg. 

The real quadratic form in the exponent of d gs , according to 60 , must be symplectically 
congruent to a diagonal one, which will turn out to be a Gibbs state. We perform the 
symplectic transformation 

(Q, P) ■= ((io/j) 1/4 q, ( W/ u)- 1/4 p); or, in shorthand, U = Su, (52) 

where S is evidently symplectic. Introducing as well the parameter A := 2-Jujfi/ (cu + /i), the 
1-quasidensity comes from the simple 

4(g,p) = -e- A(Q2+p2) . 

7T 

It is helpful to write 

A =: tanh — , whereby sinh — = - = = _ 
2' 2 v^A 2 u-n 

From the well-known series formula, valid for \t\ < 1, 



00 1 

L n {x) e~ x ' 2 t n = — e-^+t)/2(i-u 

n=0 ' 



taking t = -(1 - A)/(l + A) = -e~ p and x = 2(Q 2 + P 2 ), it follows that 



X) 



A e _ A(Q2+ p 2) = 2 ginh p y { _ iyLr{2Q 2 + 2p2) e -( Q2+ P V(2r+ i W2 _ 



7T 7T 2 

r=0 
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All formulas involving Laguerre polynomials in this section are taken from (85). We recognize 



7T 



the famous basis of orthogonal oscillator Wigner eigenfunctions on phase space 43 

f rr (Q, P) = i (-l) r L r (2Q 2 + 2P 2 ) e -^ 2+p2 \ 
Another well-known formula, 



x a L^ l (x)L c ^(x) e x dx = r -^- 5 mn , for Re a > 0, 



nl 



2.2 



guarantees the correct normalization: f f 2 r (Q,P)dQdP = (27r) _1 — see Section 

Consequently we realize that d* is indeed a thermal bath state 60 , with inverse temper- 
ature 0: 


2 



d*(q,p) = 2sinh|^e-( 2r+1 ^ 2 / rr (Q,P). 

r=0 

The occupation numbers are = 1 — e _/3 and 

n r = n e p = — v 

Clearly £ r n r = (1 - e^) £ r e~ r/3 = 1. 

We are now ready to recompute explicitly the SLK functional 

oo 

d2L(u 1 ,U 2 ) = ^ ±V n r n s frs(Ul)frs(u 2 ). 
r,s=0 

Let us write U 2 := Q 2 + P 2 . Then, for r > s, except for constant phase factors, the Wigner 



eigentransitions are known 43 



f rs {u) := - {-l) s J S -{2U 2 ) i - r - s)l2 e- i ^- s ^L r - s {2U 2 ) 



-u 2 
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T ! 



where $ := arctan(P/Q). We take f sr to be the complex conjugate of f rs . 
We sum over each sub diagonal, where r — s = I > 0: 

2^ \f n rn s frs{Ul)frs{u 2 ) 



r—s=l 



^0 1/3/2 



oo . 

e- l ^ 2 {2U 1 U 2 ) l e- a ^ + ^ e~ u "-^ ^ - L l s (2U 2 )L l s (2U 2 



7T Z 



1 



s=0 



(/ + ,)! 



7T- ' Vsinh(/3/2) / ' 

where Ii denotes the modified Bessel function, on use of 



n=0 



H] :K(x)L" n (y) = {xVt) ~ a/2 e-^M /, 



n + a)! 



1 -t 
37 



1 -t 



Similarly for r — s = — I < 0, yielding the same result replaced by its complex conjugate. 
Borrowing finally the generating function formula 



I (z)+2j2li(z) cos(Z0) = e 



z cos 6 



1=1 



where, by taking 9 = $i + $2 + 7r, one obtains for the total sum: 
1 



-(C/2+C/D/A -2l/iC/ 2 csch(/3/2)cos(i9i+i? 2 ) 



7T^ 



e -(«i+9i)(w+l»)/2 e -(p?+pi)(^- 1 + M - 1 )/2 e -gig 2 (^-M) e PiP2(^- 1 -^- 1 ) 



7T" 



which is the correct result (51), on the nose 



The choice 9 = $1 + § 2 + tt amounts to a clearly legal configuration, none other than the 
alternating sign rule for the SLK functional: = (— v/n^. Unfailingly, the analogue of 



the energy functional (46) is correctly recovered as well. Details for the latter, showing the 



correct scaling behaviour (21) in particular, are given in |86|. 



We choose to recall here that extant approximations for the exchange-correlation func- 
tional written in our terms are most often of the following form, with an obvious generaliza- 
tion of our notation: 



£xc[d] 



1 / \ f Xifc(l)Xfej(2) , 



j,k=i 



m-Q2\ 



These are all actually recipes for di- For the Miiller functional a(nj, %) = ^n^n^. A handy 



list of alternatives is provided by 87 . According to this reference, all of them (except for the 
HF functional) violate antisymmetry; nearly all of them violate the sum rule for c^; as well 
as invariance under exchange of particles and holes for the correlation part. The differences 
between Coulomb and confining potentials are of course considerable; nevertheless, analytic 
comparison of the proposed functionals with the exact one remains an useful exercise. This 
is taken up in 



Also our analysis in this section allows to throw some light from our viewpoint [89 1 on 



Gill's elusive correlation functional 16, 17, 90 



Last but not least: while our manipulations are pretty straightforward, it would be 
harder to see how to go about this problem if working on a Lagrangian hyperplane, say 
configuration or momentum space, of the full phase space. (This is somewhat analogous to 
the insight on the Jaynes-Cummings model brought about by the phase space view |91|.) 
The change of variables (52) works because there are unitary operators — the metaplectic 
representation — effecting the congruence and its inverse at the quantum level. We exhibit 



their representatives on phase space 59 , 92 



(EJ S x d* x Es)(u) = d*(Su), 



where 
2 



.(„. „, = t exp d 2 «>[(^) 1/4 -W /4 ] 



although knowledge of their existence is enough. 
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A The Thomas— Fermi workshop 



A.l Thomas— Fermi series 

Many regard the TF theory of atoms as DFT avant la lettre. Walter Kohn himself, in his 



Nobel lecture 93 , dubbed his minimum principle the formal exactification of the TF model. 



Also March has written that the "forerunner" Thomas-Fermi model is completed by the 



Hohenberg-Kohn theorem 94 . There is some exaggeration in this. But the assertions do 
apply to our minimum principle; that is to say, TF theory can justifiably be considered a 
phase space density functional theory of the type considered in this paper. 

Let us summarize the present-day model. Combining the Pauli and uncertainty prin- 
ciples, one sees that the "radius" of the atom goes like iV -1 / 3 . For a neutral atom, this 
together with ^ tell us that ground state energy should be oc —Z 7 ^ 3 hartrees. The TF 
model then gives the value cj = 0.76874512 . . . for the proportionality constant, which in 



the seventies Lieb and Simon proved (within our approximations) exact as Z f oo 95 , 96 
Nevertheless, traditional TF theory is often dismissed as accounting badly for the energy, 
being consistently larger in absolute value. Already in 1930 Dirac introduced the exchange 



correction to the TF energy 97 14 Although on the mark, Dirac's correction was not much 
used, for the good reason that it comes with a minus sign, contributing to an ever lower 
bound for the energy. 

Consider quantum mechanics (of the garden or of the Wigner-Moyal variety) for an ion 
with non-interacting electrons. The electrons would just pile up on the spectrum of states 
of a hydrogenoid system, with energies —Z 2 /2n 2 , for n the principal quantum number. 
Regarding only closed-shell ions for simplicity, and since such shells hold 2n 2 electrons, we 
exactly have 

A „ k(k + l)(2k + l) , , 9 

iV = > 2n 2 = K JK — —-t and E = -kZ 2 , 

3 

n=l 

for k closed shells. Inverting the relation gives the following rapidly converging series 

iW = (-(3/2)V3^ Ar l/3 + 1_ Z 2 _ (3^! z2Ar - 1 /3 ± ...y n 

~ (-1.145 Z 2 N 1/3 + \Z 2 - 0.073 Z 2 N~ 1/3 ± • • • ) au; (53) 
The TF variational problem without Coulomb repulsion among electrons is an exercise 



in 



iOOl Sect. 4.1], with the result 

E TF>ni = -(3/2) 1 / 3 Z 2 iV 1 /3. 

Thus the root of the problem is laid bare: the TF functional is too rough in that it "counts" 
an infinite number of states, forbidding to see anything but the leading term in the energy. 
The effect of electron screening in such a context is to reduce (3/2) 1//3 to c-j. Since the next 
term is independent of the number of electrons, it ought to be the same in the screened 



and unscreened theory, and that turns out to be the case: this is the Scott correction 101 



It is wryly amusing to note that in this paper he wrote down, prior to Wigner, the phase space quasiprob- 
ability associated to the density operator. The story has been recalled in [98]: it makes even more unsettling 



his reticence to accept Moyal's formulation of quantum theory on phase space — consult the biography 99 . 
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attributed to inaccurate treatment of the innermost shell, that was put on firmer ground 
by Schwinger |102| . Add Dirac's exchange correction and a correction term of the same 
form, worth 2/9 of Dirac's, related to the bulk electron kinetic energy and established by 

The outcome for the neutral case is the formula 



Schwinger as well 1 103 



E 



TF 



-c 7 



Z V3 + C&Z 2 _ ^^5/3 ± . . . 



au, 



(54) 



where C7 was given above, c 6 is | and c 5 ~ 0.2699. Significantly, Schwinger employed phase 
space arguments for all three terms. (For the last term, his derivation can be simplified by 
use of the native semiclassical approximation 
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Phase space approaches encapsulate 
some of the non- locality which is the bane of density gradient approximations in DFT. A 



similar observation underlies the fruitful method in 1105 



Not so long ago, the e xpression (54) was rigorously proved to be exact to the indicated 
order as Z f 00 106 107 . Due to the shell structure, no continuation at Z 4//3 order exists. 
The correlation energy is proportional to Z — see (41) and look up 94 — and so it plays 



little practical role. Empirically, the result (54) falls typically within less than 0.1% from the 



best Hartree-Fock values for ground-state energies. The model is reliable for many kinds of 
calculations, from diamagnetic susceptibilities to fission barriers. For a good review of TF 
theory, consult 
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A. 2 Making sense in WDFT of the TF scheme 

Within the old model, the energy of an atom or ion as a functional of p is given by 

E[p}=V cxt [p} + V ee [p] + T[p], 

where 



V ext [p] = ~Z J ^-dq; V ec [p] = \ J 



2J \q-q\ 



T[p] = C F I p{q) 5/3 dq, with C F 



3^4/3 
10 



and the constraint A" 



Instead one can propose the functional on phase space to be 

E[d] = V cxt [d} + V cc [d}+T[d], 

where 



(55) 
(56) 

(57) 



V exk [d] =-Z f dq, Vjd] = - ff dqdq", as before in p}; 

J \q\ 2 J J \q-q\ 

T[d] = - J \p\ 2 d(q;p) dq dp, with the constraint N = J p{q)dq. 



Here p is the known functional of d. Let us now compare the functionals in (55) with their 



partners on phase space. The TF model in the Wigner framework simply corresponds to 
the purely classical choice for A[d\. This is maybe a poor man's substitute for A[d]; but 
admissible pour les besoins de la cause. We take issue with T[p], obtained by semiclassical 
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approximation — see for instance 79, Chap. 4] or 83, Chap. 10]. This is the Schwerpunkt of 
the Thomas-Fermi approach: indeed the main source of error in the TF formalism is known 
to lie in the kinetic energy functional. In WDFT we possess instead an exact kinetic energy 
functional. In the light of the work by Dahl and Springborg on atomic Wigner functions 



reviewed in subsection |6.2[ it is clear why the uncorrected TF model fails to describe the 
strongly bound electrons, and to a lesser degree the bulk electronic gas. 

Precisely the minimum of T[p] is the value taken by T[d] if in the WDFT framework one 
settles on the ground-state phase space density 



^tf(?; p) 



Atc 3 



1 

47T 3 



e 



2Z 



x(Z 1 / 3 r/r ) - \p\ 



(5* 



d 2 x(x) x(x) 3/2 



Here is the Heaviside function, \ the solution of the famous TF equation , 

(XX \/ X 

with the boundary condition x(0) = 1, and r$ = (128/97T 2 ) -1 / 3 ~ 0.88534. By integration 
this entails p(q) = pp{q) / / 37r 2 . For the energy, invoking the virial theorem, from (58) by 
integration by parts: 



E 



e 



/2Z 
—x(Z 1 / 3 r/r ) - \p\ 



? 2 #^=fx'(o)^2 7/3 , 



It is numerically known that x'(0) — —1-5881 for the solution we are interested in. Thus the 
constant before the factor Z 7 ^ 3 is the c-j reported above. 



The question is: what error does (58) introduce? The essential point is that in the exact 



theory we are allowed to do variations only in the narrow domain of Wigner quasiprobabili- 
ties, whereas the TF equation is obtained performing an unconstrained variation (except for 
J p(q) dqdp = N), not even positivity of p being a priori necessary. This explains why the 
TF values for the energy are grossly lower than the true values. For the hydrogen ground- 
state energy, the uncorrected TF Ansatz gives Eq ~ — 0.77au, way under the — 0.5au mark. 
Since the Coulomb repulsion (that here should be put to zero) pulls the other way, this is 
more than entirely ascribable to the use of too large a set of phase space distributions, and 
to o?tf being very far from a quantum state, 



15 



Said in another way, d<?F a lousy quasiprobability makes. We remark that the associated 
electron density diverges at the origin as r _3//2 , whereas a true Wigner function must be 
everywhere continuous. (Also the behaviour of g?tf for large r differs from the oscillatory 
exponential falloff of atomic Wigner functions |62|.) In summary, the functional is good, the 



domain is bad. Or, as Parr and Yang 23 put it: "one should try to retain the good property 
of the energy functional. . . but to improve on the density" 



A. 3 The formal density matrix 



Let us, like in 21 and 



98 , formally determine the "reduced density matrix" corresponding 



to (58), by means of (9|. Introducing the relative position vector s — q — q and polar 

15 Recall that even approximating that ground state by the very small set of s-Gaussian Wigner distribu- 
tions, one obtains ~— 0.42au, now of course over the mark, a pretty better shot than the TF result. 

16 Underlying this there is the motto that functional approximation is less problematic than representa- 
bility [1091 ; this is why hope always springs eternal in generalized DFT. Our method is in this spirit. 
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coordinates for p relative to s: 



7TF(g;g 



2vr 2 



p dp l < 

o 



vr 2 s 3 .]q 



sp F (R) 



x sin x dx 



'■~ cos6 sm6d9 
1 , . 



7T 2 S 3 



SIM — X COS X j 



sp F {R) 




sin(spF(R)) spf(R) cos(spf(R)) 



ir 2 s 3 



7T 2 S 3 



MR) 



sm(spF{R)) — spp{R) cos(spf{R)) 
(spf(R)) 3 



(59) 



where we have put R := (q + q)/2. As a check, in view of 

3 5 3 5 3 5 

n / \ . 37 *C Xj 37 i-C OC 

fix) : = sin x — x cos x ~ x 1 x -\ h • • • ~ h 

JK ' 6 120 2 24 3 30 



as x 

4-0, 



we recover p(g) = 7tf(?;?)- The expression is of course ubiquitous in the theory of the 
noninteracting homogeneous electron gas — see (23| Chap. 6] or |110| ; only its manner of 
derivation is somewhat novel. 

By construction, -Di,tf is correctly normalized and hermitian; however it cannot be 
the kernel of a density matrix and it must have negative eigenvalues. As Z f oo, the 
previous expression oscillates more violently outside the vicinity of the diagonal, and the 
lower eigenvalues of 7tf migrate to zero. That means that its Wigner counterpart d^p(q;p) 
tends to become an acceptable quasiprobability; together with the relative vanishing of the 



quantum correlations 95 , this would describe how TF theory is the limit of the exact theory 
as Z f oo. Better still, one should be able to prove that g?tf approaches in some sense a 
true Wigner quasidensity (it would be most interesting to see which kind of state it looks 
like). We still ought to sharpen these remarks into a new proof of the Lieb-Simon theorem. 
Unfortunately the appropriate procedure to parametrize Wigner quasiprobability functions 
in the limit Z f oo still eludes us; but it is difficult to devise a conceptually simpler argument. 

It also stands to reason that the Scott correction and Schwinger's correction to the Dirac 
term are to be accounted for by use of the proper Wigner quasidensities. The Dirac term 



itself is easily argued for within our approach: one can read it off from (37) and (59). We 



omit this, since an identical computation is found in 23, Chap. 6]. 



A. 4 On the energy densities 



We have seen in subsection 7.5 that — as suggested already in 13 — the natural definition 



for the density of kinetic energy within phase space quantum mechanics is given by: 

1 



t 



Al\ 



(V 2 -V 2 /) 7(x;x) 



(60) 



x=x =q 



For the kinetic energy of a pure state given by a wavefunction ■?/>(<?), so that 7 = this 
leads to: 

^(g) = -^V 2 |^(g)| 2 + ^|V^(g)| 2 , 
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rather than the standard ||V?/>(<?)| 2 - Remarkably, for the ground state of a hydrogenoid ion 
this leads to a local form of the virial theorem p4] . On using (32) it is quickly seen that 



Z 4 e~ 2Zl 



2r 



The above points to good local properties of quantities obtained through the phase space 
formalism: a long standing tenet of the work by Dahl and Springborg. For instance, it 
suggests a description of the exchange hole close to being optimally localized 111 . Also for 



many atoms apparently tu obeys a local version of the Lieb-Oxford identities 112 . Now 



IVl + VH + V^V,; V 2 , 



ly 2 , 4- vH — V -V- ■ 



V 2 .. 

v r 



Then with x = Pf(R)s and with the help of (59) and (60) we compute 



t M {q) 



3 -p(q)(d 2 s + - s d s )f(x) 



C F p(q) 



5/3 



s=0 



We have checked that T[<i TF ] = T[p TF ]. However, we are walking on a tightrope. For 



reconsider (60), under the more general form 



a 



x=x =q 



where a + (5 = ^, in the framework of spectral expansions; it yields 



(Vp.) 2 g 
8p* 2 



V 2 A 



El 



(Vv^) 2 -/3V 2 Pi 



put a 

obtained instead t 



| holds only for true 
= |; suppose we had 



The formula with any coefficients a, (3 respecting the sum rule a + (3 = 
quantum states. The phase space formalism seems to "choose" a = (3 

0, forcing /? = |, in the calculation (60). Then we see at once that we would have 

= -|v|p(g) + c FP (qf 3 



This is all very well for a real atom, since 
then f V| p(q) dq = 0; but for the Thomas-Fermi model J V? p(q) dq diverges. 

The appearance in the last display of the term of the von Weizsacker type, everywhere 
positive and finite, is welcome. This kind of term is known to tame the worse aspects of 
the TF functional, in particular allowing for a good behaviour at the atomic nucleus site j^] 
Because of the Schwarz inequality, the von Weizsacker functional W[p] is convex; this means 
in particular that 

(VP) 2 f(V Pi 



W[p] :- 



dq, 



17 The book by Parr and Yang contains an interesting development in this connection: it derives the 
improvement of the TF energy functional by such terms from the semiclassical expansion of the Wigncr 
function [23j Sect. 6.7]. However, the coefficients obtained for the gradient corrections are not optimal, and 
such expansions, when pushed to higher order, open their own box of horrors in the form of hopelessly 
divergent contributions. 



43 



for pi as in the foregoing; for an atomic system the integral on right hand side should precisely 
give the total kinetic energy 



IN 



It is interesting to revisit harmonium in this context. It is not hard to see 52 that the 
kinetic energy density for this system is of the form 



t(q) = W[p(q)]+3p(q) 



(u - pf 



S(u + p) 

Note that only the von Weizsacker term contributes in the weak repulsion limit. We run a 
check on this. Enter equation (26); from it we glean that 

2 



W[p sa ](q) 



up 
CO + p 



Q 2 pgs(q)- 



Therefore 



T = / W[ Pgs \(q) d A q 



,, | 3 (a; — p) 2 3up | 3 (a; — p) 2 3(u + p) 



4(cu + p) uj + p A{uj + p) 



as it should be. 



B Characterization problem for Wigner functions 



Nothing prevents making variational calculations with the functional (57), or appropriate 
corrections of it, by trial families of Wigner Gaussian orbitals, using the formulas of subsec- 
tion 16.11 The facts indicated at the end of subsection IA.1I ensure reasonable results. 

Those who do not fancy Gaussian basis sets then face the question of characterizing 
Wigner quasiprobabilities. A one-body Wigner function is just N times a convex combination 
of pure states on a single copy of phase space. Thus the question is how to recognize a 
quantum pure state representative among all functions on phase space. Beyond reality and 
square summability, some necessary conditions are easy: continuity (not smoothness, as is 
sometimes erroneously assumed); the bound ([6]); positivity of the integrals ([8]), and indeed of 
the integral over any Lagrangian plane of the phase space. A less known necessary condition 
is the following. Given a wave function for each q, p consider the translated wave function 



ip -X 



ty{x — q) e 

Its Wigner counterpart is clearly given by P^(q — q',p — p), if P-q, corresponds to Thus 
we also have the necessary condition 



/ 



M(q; p)P(q — q;p — p ) dqd 3 p > 0, 



for all pure states P, all Wigner quasiprobabilities M and all q,p. 

Necessary and sufficient conditions are harder to come by. The defining condition 
P = (2tt) 3N (P x P) is difficult to handle. Equivalent conditions will have an oscillatory 



18 



Convexity of W links it to relative entropy; but we refrain from going into that. 
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integral somewhere; however, the following alternative may appear more convenient to some. 
Performing a Fourier transform on the second set of variables, we obtain 

P(q-r):= J P(q- p)e~ 2i ^ dp = J J V(q - z)V(q + z)e 2i ^~ r) dz dp 

Thus, with P(s) := P(s; s), we get 

P ( | (9 + ?)) = *(0)*(« + f) I P(!(9-?)) = pip*(<z-r)*(0); 
P(0) = ^*(0)*(0); » that Mi±|MiZl)) = 

This condition is easily seen to be sufficient as well. Also note that 

A g log P{q;r) = A r log P{q;r). 
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